Symbolic Techniques in MathModelica System Designer

How to achieve efficient simulation of Hybrid DAEs  from Modelica models.

Peter Aronsson
Product Manager

MathCore Engineering AB, Linköping Sweden

Abstract

MathModelica System Designer is a modeling and simulation environment that allows creation of mathematical models in an intuitive and easy manner by using the object oriented, equation based language Modelica. This is achieved in a graphical model editor where components are connected together to form models of systems. Since Modelica is equation based it is suitable for any domain, and it is therefore also possible to mix components of different domains in the same model, like for instance mechanical, electrical or discrete controller components.
A Modelica translator is used to translate a Modelica model description into a set of differential and algebraic equations with discrete events, a so called hybrid DAE system. Due to the object oriented nature of the Modelica models such DAE systems are typically sparse and often have high index. The index of a DAE corresponds to how many times certain equations must be differentiated in order to transform the DAE into explicit ODE form. In order to achieve efficient simulation of such models it is often required to symbolically simplify the equation system and reduce its size. Also, to be able to use an arbitrary numerical solver for the model it is also required to reduce the index of the system by a technique called index reduction.

In this presentation we show how these techniques are used in MathModelica System Designer to achieve efficient simulation of Hybrid DAEs. We also present how to utilize this to be able to solve Hybrid DAE systems in Mathematica by using the Mathematica connection in MathModelica System Designer.

Introduction

A Modelica model is built in MathModelica System Designer by drag and drop of components and connecting these togheter to form a system.

Graphical modelling

MathModelica System Designer has a graphical model editor where models can be built using drag and drop in an easy manner. A DCmotor consisting of an electrical part and a mechanical part can be produced by using components from the Modelica standard library.

The figure above shows the MathModelica Model Editor and the DCMotor example. It consists of an electrical part to the left and a mechanical part to the right. The different components has a strong correspondence to real physical objects, which is not the case for a  block-based model of the same system.

The above model can be simulated in Mathematica by using Simulate.

<< MathModelica`

Simulate[GettingStarted . MultiDomain . DCMotor, {t, 0, 10}]

This translates the model into simulation code and executes it. The result of the simulation is then loaded into Mathematica and can be plotted using PlotSimulation.

PlotSimulation[{inertia1 . w[t]}, {t, 0. , 10. }];

Translating to Hybrid DAE

The first step in translating the Modelica model into a runnable simulation is to produce the complete equation system of the models. This is performed by a process called instantiation. All the components of the models are instantiated and the variables and equations are collected to form the complete system. Also, equations from the connections between components are formed. The result is a list of n variables and a list of n equations.
The textual description of these variables and  equations is presented below.

Textual description of the equation system of a model

SimulatorCommand["instantiateModel(GettingStarted.MultiDomain.DCMotor)\
"]

Quit[];

In this section we present the different techniques applied to the equation system in order to produce an efficient simulation.

Elimination of simple equations

Since the equation system from Modelica models is typically sparse due to for instance the equations generated from connection of different components. Thus, much is gained by elimination of these variables. The translator will remove equations on the form

a  ==  b,  c  ==  -d,  e  =  4.5

This will typically reduce the system by several magnitudes. This will also speed up later steps of the translation process.

Transforming to explicit ODE

The next step of the translation is to transform the DAE into an explicit ODE form, at least conceptually. In order to do this the equations are sorted into a so called BLT (Block Lower Triangular) form. The goal of this process is to produce a function

\(TraditionalForm\`X\& .  = f(X, Y, t)\)

where X are the state variables and Y are the algebraic variables and t is the independent time variable. This function can contain systems of equations that needs to be solved simultaneously. Such systems can be both linear (which can in some cases be solved during translation by gaussian elimination) or non-linear.

The advantage of transforming the equation system to explicit ODE form is that once this is done, any numerical solver can be used, including for instance explicit runge kutta solvers that only works for ODEs. This is particularly useful in real time simulations where fast, fixed step-size solvers are required.

Index reduction using dummy derivatives

Definition of Index

The Pantilides algorithm

For high index problems it is necessary to reduce the index to be able to transform the DAE to explicit ODE form. A high index problem occurs when there are constraint between state variables.
As an example, let us consider two connected inertias.

Off[General::"spell1"];

The equations of the model above are:

eqns={(Subscript[ω, 1])[t]==((Subscript[ϕ, 1]))'[t],(Subscript[a, 1])[t]==((Subscript[ω, 1]))'[t],(Subscript[ω, 2])[t]==((Subscript[ϕ, 2]))'[t],(Subscript[a, 2])[t]==((Subscript[ω, 2]))'[t],(Subscript[J, 1]) (Subscript[a, 1])[t] + (Subscript[J, 2]) (Subscript[a, 2])[t] == k,(Subscript[ϕ, 1])[t] == (Subscript[ϕ, 2])[t]}

params  =  {(Subscript[J, 1]) -> 0.1, (Subscript[J, 2]) -> 0.2, 
          k -> 2};

sol = 
      NDSolve[Join[
          eqns /. params, {(Subscript[ω, 1])[0] == 0, (Subscript[ω, 2])[0] == 
              0, (Subscript[ϕ, 1])[0] == 0, (Subscript[ϕ, 2])[0] == 
              0}],  {(Subscript[ω, 1]), (Subscript[ω, 2]), (Subscript[ϕ, 1]), (Subscript[ϕ, 2]), (Subscript[a, 1]), 
          (Subscript[a, 2])},  {t,  0,  1}]

NDSolve[
      Join[eqns /. 
          params, {(Subscript[ω, 1])[0] == 0, (Subscript[ω, 2])[0] == 
            0, (Subscript[ϕ, 1])[0] == 0, (Subscript[ϕ, 2])[0] == 
            0}],  {(Subscript[ω, 1]), (Subscript[ω, 2]), (Subscript[ϕ, 1]), (Subscript[ϕ, 2]), (Subscript[a, 1]), 
        (Subscript[a, 2])},  {t,  0,  1}, Method -> "BDF"]

Plot[(Subscript[ω, 1])[t]  /. First[sol], {t, 0, 1}];

Lets give a more interesting input:

params2  =  {(Subscript[J, 1]) -> 0.1, (Subscript[J, 2]) -> 0.2, 
        k -> Sin[2 π t]}

sol2 = 
      NDSolve[Join[
          eqns /. params2, {(Subscript[ω, 1])[0] == 
              0, (Subscript[ω, 2])[0] == 0, (Subscript[ϕ, 1])[0] == 
              0, (Subscript[ϕ, 2])[0] == 
              0}],  {(Subscript[ω, 1]), (Subscript[ω, 2]), (Subscript[ϕ, 1]), (Subscript[ϕ, 2]), (Subscript[a, 1]), 
          (Subscript[a, 2])},  {t,  0,  10}]

Plot[
        Evaluate[{(Subscript[ω, 1])[t], (Subscript[ω, 2])[t]}  /. First[sol2]], {t, 0, 
          10}];

Plot[
      Evaluate[(Subscript[ϕ, 1])[t] - (Subscript[ϕ, 2])[t]  /. First[sol2]], {t, 0, 10}, 
      PlotRange -> All]

Using pantilides to reduce the index

The pantilides algorithm reduces the index of the dae by differentiating equations. The inertias system can be reduced by differentiating the constraint equation (Subscript[ϕ, 1])[t] == (Subscript[ϕ, 2])[t]:

(eqn1 = D[eqns[[-1]], t])

By introducing the new equation to the system we have reduced the index by one.

This system can now be solved for the first order derivatives. We form the new system by replacing the constraint equation by the differentiated one.

newEqns = Join[eqns[[{1, 2, 3, 4, 5}]], {eqn1}]

We eliminate (Subscript[a, 1]) and (Subscript[a, 2]): ( why doesn't Eliminate work here?)

newEqns2  =  
      newEqns  /.  {Rule  @@ newEqns[[2]],  
            Rule  @@ newEqns[[4]]} // DeleteCases[#, True] &

Equation 1,2 and 4 still give a constraint between states, differentiate these and form a new system:

newEqns3  =  {newEqns2[[3]], D[newEqns2[[1]], t], 
        D[newEqns2[[2]], t], D[newEqns2[[4]], t]}

This can now be solved using an ODE solver, i.e. the index of the system has been reduced to zero.

\(sol3 = 
      NDSolve[Join[
          newEqns3 /. 
            params2, {ω\_1[0] == 0, ω\_2[0] == 
              0, ϕ\_1[0] == 0, ϕ\_2[0] == 
              0, \(ϕ\_1'\)[0] == 0, \(ϕ\_2'\)[0] == 
              0}], \ {ω\_1, ω\_2, ϕ\_1, ϕ\_2}, \ {t, \ 0, \ 200}, Method -> "BDF", MaxSteps -> 50000]\)

Plot[
        Evaluate[{(Subscript[ω, 1])[t], (Subscript[ω, 2])[t]}  /. First[sol3]], {t, 0, 
          10}];

However, the constraint between is not enforced by the solver, thus causing the solution to drift:

Plot[
        Evaluate[(Subscript[ϕ, 1])[t] - (Subscript[ϕ, 2])[t]  /. First[sol3]], {t, 0, 
          200}];

The dummy derivatives method

To avoid the drift problem when reducing the index a method called the dummy derivates method can be used:

When the pantilides algorithm detects an higher index problem, the constraint equation(s) is/are differentiated and added to the system. One of the derivatives is transformed to an algebraic variable, a so called dummy derivative.

We again look at our Connected inertias example. The first equation to become differentiated is:

(eqn1 = D[eqns[[-1]], t])

The dummy derivatives method selects one of these variables as a state and transforms the other to an algebraic variable. We select (Subscript[ϕ, 1]) as state variable and transform ((Subscript[ϕ, 2]))' to a dummy derivative (and (Subscript[ϕ, 1]) then becomes a dummy state), resulting in the new system:

newEqns4 =Join[eqns, {eqn1}] /. {((Subscript[ϕ, 2]))'->d(Subscript[ϕ, 2])}

Again,  we eliminate (Subscript[a, 1]) and (Subscript[a, 2]) for convenience

newEqns5  =  
      newEqns4  /.  {Rule  @@ newEqns[[2]],  
            Rule  @@ newEqns[[4]]} // DeleteCases[#, True] &

NOTE: The equation (Subscript[ϕ, 1])[t] == (Subscript[ϕ, 2])[t]is no longer a constraint between two states.

The next step will differentiate the three other equations giving constraint between states ((Subscript[ω, 1])[t]==((Subscript[ϕ, 1]))'[t],(Subscript[ω, 2])[t] == d(Subscript[ϕ, 2])[t],((Subscript[ϕ, 1]))'[t]==d(Subscript[ϕ, 2])[t]) as previously:

newEqns6 =  
      Join[newEqns5, {D[newEqns5[[1]], t], D[newEqns5[[2]], t], 
          D[newEqns5[[5]], t]}]

We select (Subscript[ω, 2]) as state and transform ((Subscript[ω, 1]))' to a dummy derivative (and other derivatives to algebraic vars)

newEqns7=newEqns6 /. {((Subscript[ω, 1]))'->d(Subscript[ω, 1]),(d(Subscript[ϕ, 2]))'->dd(Subscript[ϕ, 2]),((Subscript[ϕ, 1]))''->dd(Subscript[ϕ, 1])}

sol4 = 
      NDSolve[Join[
          newEqns7 /. 
            params2, {(Subscript[ω, 2])[0] == 0, (Subscript[ϕ, 1])[0] == 
              0}],  {(Subscript[ω, 1]), (Subscript[ω, 2]), (Subscript[ϕ, 1]), (Subscript[ϕ, 2]), 
          d(Subscript[ϕ, 2]), d(Subscript[ω, 1]), dd(Subscript[ϕ, 2]), dd(Subscript[ϕ, 1])},  {t,  0,  
          200}, MaxSteps -> 50000]

Plot[
        Evaluate[{(Subscript[ω, 1])[t], (Subscript[ω, 2])[t]}  /. First[sol4]], {t, 0, 
          10}];

We now see that the drift phenomenon has been eliminated:

Plot[
        Evaluate[(Subscript[ϕ, 1])[t] - (Subscript[ϕ, 2])[t]  /. First[sol4]], {t, 0, 
          20}];

In comparison with using only pantilides:

Plot[
        Evaluate[{(Subscript[ϕ, 1])[t] - (Subscript[ϕ, 2])[t]  /. 
              First[sol3], (Subscript[ϕ, 1])[t] - (Subscript[ϕ, 2])[t]  /. 
              First[sol4]}], {t, 0, 200}, 
        PlotRange -> All];

Quit[];

Example in MathModelica System Designer

<< MathModelica`;

MathModelica System Designer performs index reduction using dummy derivatives automatically. This is necessary because:

The Model used in the example above can easily be made using the Model editor.
Start Model Editor

PrintDefinition[Model1]

Simulate[Model1, {t, 0, 200}, 
      NumberOfIntervals -> 50000]

PlotSimulation[inertia1 . w[t], {t, 0, 20}];

PlotSimulation[
        inertia1 . ϕ[t] - inertia2 . ϕ[t], {t, 0, 
          200}];

Summary

Download at:  http://www.wolfram.com/products/applications/mathmodelica